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Formation of a sonic horizon in isotropically expanding Bose-Einstein condensates 
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We propose a simple experiment to create a sonic horizon in isotropically trapped cold atoms 
within currently available experimental techniques. Numerical simulation of the Gross-Pitaevskii 
equation shows that the sonic horizon should appear by making the condensate expand. The ex- 
pansion is triggered by changing the interaction which can be controlled by the Feshbach resonance 
in real experiments. The sonic horizon is shown to be quasi-static for sufficiently strong interaction 
or large number of atoms. The characteristic temperature that is associated with particle emission 
from the horizon, which corresponds to the Hawking temperature in an ideal situation, is estimated 
to be a few nK. 

PACS numbers: 03.75.Kk, 03.75.Hh, 04.62.+V, 05.30.Jp 



I. INTRODUCTION 



For the exploration of cosmology and gravitational 
physics, it is necessary to have a deep understanding of 
quantum filed theory in curved spacetime: It is widely 
believed that everything except for the spacetime it- 
self should originate from quantum fluctuations in the 
early Universe. Quantum effects on curved spacetime, 
such as the Hawking radiation, give us theoretical sup- 
port for black hole thermodynamics. However, it is ex- 
tremely hard to verify such quantum effects experimen- 
tally. For instance, the Hawking radiation is thermal 
radiation emitted from a dynamically formed stationary 
black hole However, the characteristic temperature 
of the thermal radiation, the Hawking temperature, is on 
the order of several tens of nanokelvins at most, which is 
much lower than the cosmic microwave background radi- 
ation temperature. So detecting thermal radiation from 
a real black hole is almost impossible. 

One way to circumvent this difficulty is to make use of 
artificial black holes [|J • Unruh showed in his seminal 
paper Q that excitations in a supersonic flow corresponds 
to a scalar field equation on a curved spacetime includ- 
ing a horizon. Since the phenomenon of the Hawking 
radiation can be separated from gravitational physics, it 
is possible to detect the corresponding phenomenon in a 
fluid system with sonic horizon The basic idea is to 
identify fluid flow with curved spacetime and excitation 
modes with fields on the curved spacetime. A black hole 
event horizon corresponds to a sonic horizon in a fluid. 
For the purpose of investigating the quantum effects, a 
quantum fluid should be considered. As such a quantum 
fluid, Bose-Einstein condensates (BEC) in trapped cold 
atoms @,@| are one of the most suitable systems 
A crucial advantage is that one can control scattering 
length between atoms by making use of the Feshbach 
resonance [10|. In fact, that experimental technique was 
used in observing jets and bursts in a collapsing conden- 
sate, which is called "Bose- Novae" [ll|. An remarkable 



explanation of burst and jet phenomena in Bose-Novae 
was proposed in 12j lSif. based on quantum field theory 



of particle creation and structure formation in cosmolog- 
ical spacetime. 

In order to verify the Hawking effect in fluid analogy, it 
is necessary to create a stationary sonic horizon because 
it is a phenomenon on a dynamically formed stationary 
black hole. Although several possibilities have been dis- 
cussed so far @, [1] , it seems difficult to realize exactly sta- 
tionary sonic horizon in cold atoms. However, if one can 
make a quasi-static horizon for high frequency modes, 
particle emission from the horizon is also expected. In 
this paper, we numerically demonstrate that a quasi- 
static horizon is realized without introducing new exper- 
imental techniques beyond currently available ones. We 
consider an expanding BEC driven by a sudden change of 
the interaction. Numerically solving the Gross-Pitaevskii 
(GP) equation, we show that a quasi-static horizon ap- 
pear. 

We note that there have been made great efforts to cre- 
ate cosmological geometry using expanding BEC 14|-[2l|. 
In these papers, the analogue models with specific cos- 
mological metrics such as Friedmann-Robertson- Walker 
(FRW) metric or de Sitter metric were discussed and the 
effects of particle creation in these cosmological space- 
times were investigated. However, in this paper, we do 
not intend to obtain any cosmological analogue model 
with well-known analytic metric. But we try to obtain 
dynamically formed quasi-static sonic horizon. Further- 
more, the sonic horizon should be formed in hydrody- 
namic regime of the condensate because the spacetime 
analogy is only valid in such regime. The appearance 
of a horizon due to expansion of a condensate was no- 
ticed in the previous works, for example in [l5l |. and its 
formation itself is not surprising. But it is non-trivial 
whether the condensate flow at the horizon is in the hy- 
drodynamic regime, or not. In this paper, we show that 
the quasi-static sonic horizon will appear in the hydrody- 
namic regime of the condensate by changing the atomic 
interaction instantaneously. 
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II. ANALOGUE SPACETIME IN BEC 

In the coherent state path integral formulation, the 
action of bosons is given by 
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--UolW), (1) 

where V ex ^ is the confining potential and Uq = 4irh 2 a/m 
with a the s-wave scattering length. For (f> and V ex t, 
spatial and time dependences are implicit. The saddle 
point equation for this action leads to the GP equation: 



m t * = ( -^v 2 + U cxt + C/ |*| 2 ) 4- 



(2) 



This GP equation governs the dynamics of the conden- 
sate whose order parameter is given by ^ . 

Now we consider hydrodynamical approximation. We 
denote the bosonic field <\> as <j> — VPo + pe l ^ vo+ ' p \ where 
y/po and tpo are the amplitude and the phase of re- 
spectively. (Namely, 4 = ^/po exp(«^ )-) The fields p 
and (p describe the non-condensate part of the bosonic 
field. If the density gradient is sufficiently smooth over 
the scale determined by the local healing length £(r, t) = 
h/(2nipoUo) 1 ^ 2 , or, in other words if the conditions, 



|£Vp /po| 2 «l and |£Vp/p| 2 «l, 



(3) 



are satisfied, hydrodynamical approximation is justified. 
(The condition shall be examined later.) Under the 
above condition, the equation for p is 



p = -h(<p + v • Vtp)/U , 



(4) 



where Vo = -^Vfo is the background fluid velocity, and 
the effective action for ip is 
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(ip + v • Vtp) 2 



^(V^) 2 



.(5) 



Taking variation with respect to ip, we find that the field 
equation for tp has the form of a propagating wave equa- 
tion. Also, the equation for the field ip can be expressed 
as d/j.(V~99 IJ ' 1 ' 'd„)ip = 0, where g^ v is the inverse of the 
following matrix: 



9^v oc 



v 



-v 

1 



(6) 



with c s = \J poUo/m and g = detg^. Thus, the equa- 
tion is equivalent to an equation for a massless field on 
a curved spacetime determined by the metric (JSJ with c s 
the speed of "light." Note that in order to interpret the 
quantity c s as a velocity, Uq must be positive because, for 
negative Uq, c s becomes pure imaginary. Hereafter we 
consider positive Uq, which leads to an effective space- 
time with Lorentzian signature. 



For the excitation modes of tp whose frequencies, say 
u>, are much higher than the frequency ujbec, which is 
associated with the condensate motion, the condensate 
will be quasi-static. (For moderate changes of the inter- 
action, wbec turns out to be the trapping harmonic po- 
tential frequency u>hoi as shall be discussed below .) The 
analogy between fields on the curved spacetime and ex- 
citation modes on the fluid flow is meaningful only when 
the conditions ([3]) are satisfied. The latter condition in 
Eq.© turns out to be u? < (c s /£) 2 , by using Eq.®. 
Thus, the frequency uj has an upper limit. The former 
condition in Eq.Q is satisfied in the regions far from 
the edge of the condensate. (In contrast, if one is very 
close to the edge, zero-point oscillations become domi- 
nant, and so the former condition in ([3]) is not satisfied.) 
If there exists intermediate region for ui of 



wbec < w < c s /£, 



(7) 



then the hydrodynamical approximation is justified and 
the condensate is quasi-static for excitation modes. Note 
that those modes are associated with particle emission 
from the horizon if the hydrodynamical flow has a dy- 
namically formed sonic horizon. The necessary condition 
for the existence of the intermediate region ([7]) is 



> 1. 



(8) 



In the following, we mainly consider condensate satisfy- 
ing the above condition. 



III. FORMATION OF SONIC HORIZON 

Now we investigate sonic horizon formation in an ex- 
panding BEC trapped in isotropic harmonic potential, 
V cx t = m < jJ \ r 1 j2, where r is the radial coordinate. Ini- 
tially, we set the condensate in a ground state with an 
initial atomic interaction a^. At t — 0, the atomic inter- 
action is changed suddenly from a, to a/(> <2j), which 
makes the condensate expand. Then, formation of sonic 
horizon can be expected. The reason is as follows: The 
sound velocity is proportional to square root of the con- 
densate density and a decreasing function of r. In con- 
trast, the fluid velocity is an increasing function of r 
and the condensate expands fast around its edge whereas 
v (r = 0) = due to the boundary condition. Therefore, 
at an intermediate radius, vq exceeds c s and the fluid flow 
is transonic. It has a surface satisfying c s = \vq\ which 
is called a sonic horizon. We should note that the sonic 
horizon corresponds to a horizon in the analogue space- 
time defined by the metric ([6]). We also note that the 
subsonic region is around the center of the condensate 
and inside of the sonic horizon. 

In general, if a fluid has a static sonic horizon and a 
proper quantum state for an excitation field is realized, 
then it is theoretically predicted that the horizon will 
emit thermal radiation of the quantum field. As will 
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FIG. 1: Times t c in each simulation are shown in units of 
^ho ■ The horizontal axis is the ratio of a/ to en. 

be discussed in Appendix if the sonic horizon in the 
expanding condensate is quasi-static for the field ip, the 
horizon will emit thermal radiation into the center of the 
condensate. The temperature characterizing the thermal 
emission (Hawking temperature) is given by the following 
formula: 

t p c = Tir dr ( Vo ~ Cs )l r H> ( 9 ) 

zttkb 

where u_is the horizon radius and is the Boltzmann's 
constant |j]|24(. From the above expression, it is found 
that the Hawking temperature is determined by gradi- 
ent of fluid and sound velocity at the horizon. Thus, 
it is important to investigate the velocity gradients at 
the horizon. In the derivation of the formula © , it is as- 
sumed that the dynamically formed horizon is static, but 
in actual experiments, this assumption is not satisfied ex- 
actly. Therefore, the spectrum of the particle is not fully 
given by the single Planck's distribution function, but 
rather given by a superposition of the Planck's distribu- 
tion functions with slightly different temperatures. Even 
if this is the case, the energy scale of the particle creation 
emitted from the dynamically formed horizon is on the 
order of T pc . 

We have simulated the expansion of the condensate by 
solving numerically (using the Crank-Nicolson scheme) 
the time-dependent GP equation. The initial ground- 
state wave function is obtained by solving the GP equa- 
tion using the steepest descent method for an initial s- 
wave scattering length ai and the number of atoms N . 
We have computed c s and the radial velocity of the con- 
densate via 

c s = y/(W)U /m, (10) 
v = h[V*d r V- (a r **)*]/(2mi|*| 2 ), (11) 

and searched for parameter sets leading to |uq| > c s- 

In the following, we assume that the condensate con- 
sists of N — 10 5 Rb atoms. (The values of atomic inter- 
action given below are those in the case of N = 10 5 . If 
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FIG. 2: Sound velocity c s (solid line) and the fluid velocity 
vo (dashed line) versus r at t = 0.4 in the case of a,i — 200ao 
and a/ = 5ai are shown in units of (TioJho/m) 1 ^ 2 ■ The 
healing length ^(dotted line) is shown as well in units of 

aho = (h/muj ho ) 1/2 . 

N = 10 5 /n with an integer n, then ai and a/ should 
be multiplied by n.) The initial atomic interaction 
is assumed to be tn = 50ao, 200ao and 800ao where 
a = 0.53 x 10~ 10 m is the Bohr radius. The follow- 
ing change of the atomic interaction has been simulated: 
aj/at = 2,3,4,5,6,7,8,9,10,15,20,25,30. 

Just after t — 0, the condensate begins to expand in 
the trapping potential and the expansion is accelerated 
for a while. At some time, say t = t c , the expansion 
turns to be decelerated. Figure Q] shows t c as a function 
of a/ /dj . It is seen that t c does not depend on the initial 
strength of the interaction. For t := Uhot < tt/2, the 
condensate continues to expand, and at t ~ tt/2, the 
condensate starts to collapse. Therefore, we turn off the 
trapping potential at t = 7r/4 and make the condensate 
expand freely in order to keep the horizon for a while. 

As far as we have investigated, sonic horizon always 
appears in the sence of the surface where vq exceeds c s . 
As an example of a sonic horizon, we show Fig. [2] which 
is the snapshot at t = 0.4 in the case of a, = 200ao and 
a/ = 5a;. We see that, around r = 7a^ Q , the fluid velocity 
exceeds the sound velocity, and the sonic horizon exists 
there. In this case, we find that, at t = u>h t = 0.11, 
the horizon appears. Fig. [3] shows the time dependence 
of the radius of the horizon, say m, and the velocity 
gradient at the horizon d r (vo — c a )\ rH . 

If we keep the trapping potential for a long time, an 
oscillating behavior of the condensate is observed. The 
period of the oscillation is about n in units of wZ Q and 
we find o>bec — ^ho, within a moderate change of the 
interaction. This oscillation is just like an oscillation of a 
droplet confined in a harmonic potential. Therefore, the 
condition © can be rewritten as 

7^ » 1- (12) 

C,(jJho 
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FIG. 3: Time dependence of d r (i)o — c s )r H (solid line) in units 
of ujho and the position of the sonic horizon ru (dashed line) 
in units of auo = (fi/muiho) 1 ^ 2 in the case of a, = 200ao and 
aj = 5ai. 



FIG. 5: Horizon life time as a function of afja% is shown in 
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FIG. 4: c a /(,uiho as a function of a// at is shown for each 
initial scattering length. 



Now, we are interested in sonic horizon where (fl~2|) is sat- 
isfied. The intermediate region ([7]) exists if, for example, 
the following inequality is satisfied: 



> 22.5. 



(13) 



For this choice of the lower bound, there exists a re- 
gion of w satisfying both conditions of u> > 10w/ lo and 
(c s /w£) 2 > 5. The condition (fT3|) ensures hydrodynamic 
flow and quasi-static nature of the condensate. Fig. 0] 
shows c s /^uJho as a function of a//<Zj. We define horizon 
life time as the time interval during which the condition 
(fl~3"l) continues to be satisfied at the horizon. The horizon 
life time is shown in Fig. [5] As far as we have investi- 
gated, the condensate flow satisfying (fT3)) at sonic horizon 
appears only when af > 5a t for a, = 50a , at > 4a.j for 
a, : = 200a and a/ > 3a^ for a.j = 800a . 

The Hawking temperature at t = t c and t = 0.79 (just 



units of w h * . 
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FIG. 6: Hawking temperature at t = t c in units of nK. In 
the evaluation, we assume Uho = 1400 Hz. 



after turning off the trapping potential) are shown in Fig. 
[6] and Fig. [3 respectively. In the evaluation, we assume 
the frequency u>ho = 1400 Hz. The Hawking temperature 
at t = t c depends on the ratio a//a, almost linearly. In 
contrast, for a//aj > 9, the Hawking temperature at 
t = 0.79 does not depend on the ratio so much. From the 
simulations, the temperature is expected to be a few nK. 

For this spherically symmetric trap, one may concern 
the three-body recombination loss of condensed atoms. 
Now, we check the effect of three-body losses for the given 
peak density. This effects may be taken into account by 
incorporating the imaginary term describing the inelastic 
process in the GP equation [22j 



2m 



where K 3 denotes three-body recombination loss-rate co- 
efficient. Then, the three-body loss is proportional to the 
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FIG. 7: Hawking temperature at t = 0.79 in units of nK, 
when just after the trapping potential is turned off. In the 
evaluation, we assume luho = 1400 Hz. 




FIG. 8: Who-dependence of Th and horizon lifetime in the 
case of ai = 200ao and a/ = lOa^. The Hawking temperature 
is shown in units of nK and horizon lifetime is in units of co^ . 



cube of the atomic density 
d 



at 



m\ 2 d A r = -K 3 / |*| b rfV, 



which implies that the three-body loss rate is given by 
R 3 = K 3 J\^\ 6 d 3 r/ J \^\ 2 d 3 r. For the value of K 3 , we 
assume K 3 — 2 x 10~ 28 cm 6 /s, according to [13]. Of 
course, high atomic density causes many inelastic pro- 
cesses and gives high atomic loss rate. In our numerical 
simulation, the upper limit of the loss rate can be esti- 
mated by use of the peak density as R 3 < 3 x 10s _ 1 , where 
the total atomic number was set to be N = J |<F| 2 d 3 r = 
10 5 . Then, the three-body loss can be ignored because 
we consider the time scale of < 10 ms. 

In the above evaluation for Hawking temperature, hori- 
zon lifetime and R 3 , we have assumed that the trapping 
frequency is U} lo = 1400 Hz. Note that w/, is the en- 
ergy scale of the system. Therefore, a large value of 
u>ho is plausible to increase the characteristic tempera- 
ture for the particle emission, though the time evolution 
process becomes rapid for large loho- If lower frequency 
is assumed, lower temperature, longer horizon lifetime 
and fewer three-body loss rate would be expected. As 
an example, Fig. [8] shows w^o-dependence of the Hawk- 
ing temperature and the horizon lifetime in the case of 
a, = 200a and a,f = lOa^. 



IV. BOGOLIUBOV SPECTRUM 

In the above numerical simulations, we assume there 
is no dynamical instability. Now, we check whether there 
is dynamical instability or not, within Gaussian approx- 
imation. For that purpose, we study the Bogoliubov-de 
Gennes equations: the second quantized field equations 



for the excitation helds 6cj> and S<fi are given by 

h 2 



ihdf 



2m 
fi 2 



X7 2 + V cxt + 2U \»£\' <ty + E/b* 



ihd t 6(p = -—V 2 + V c 
2m 



2u \m 



UofW) 



The excitation spectrum is computed by performing the 
Bogoliubov transformation: 

H = J2 { u « (r) b a e- iE ^ h - v a (r) fit e ^-t/»j >(14) 



■iE a t/h 



.(15) 



The energy spectrum E a is calculated by diagonalizing 
the skew symmetric matrix, which is carried out by using 
a routine in LAPACK. For the parameter values taken 
above, we find that all eigenvalues do not have the imag- 
inary parts within numerical errors. Therefore, within 
Gaussian approximation, there is no dynamical instabil- 
ity. In addition, we find that there is no level crossing. 



V. SUMMARY 

To summarize, we have proposed an experiment to 
create a quasi-static sonic horizon using an expanding 
BEC. It has been shown that the dynamically formed 
quasi-static sonic horizon is in hydrodynamic regime as 
it should be to discuss analogy with curved spacetime in 
BEC. Under suitable choices of the interaction parameter 
and the confining potential, the characteristic tempera- 
ture of the particle emission is expected to be a few nK 
for sufficiently strong confining potential. Large num- 
ber of atoms or strong atomic interaction improves the 
quasi-static nature of the horizon. 
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Of course, other effect such as cosmological particle 
creation can occur in this expanding BEC setup, as dis- 
cussed in [3]-[2l|. In this paper, we have focused on 
how to make dynamically formed quasi-static sonic hori- 
zon in the hydrodynamic regime of the condensate flow. 
In order to investigate cosmological particle creation ef- 
fect and other excitations arising from depletion, we need 
a different numerical simulation scheme. The result will 
be reported in a future publication. 

Furthermore, it is interesting to investigate numeri- 
cally the behavior of negative frequency modes with pos- 
itive norm which seem to be related to Hawking effect as 
was discussed in 2^] [27] . This point shall be investigated 
in a future publication. 
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APPENDIX A: PARTICLE CREATION 
PHENOMENON 

Here we focus on spherically symmetric quantum fluc- 
tuations by symmetry. At t < 0, the fluid velocity vq = 0, 
and the metric of the initial static effective spacetime is 



ds 2 oc —c 2 dt 2 



dr 2 



l dni 



(Al) 



where dVl 2 s2 is the element of solid angle on the unit 
sphere S 2 . After the increase of the interaction, the ef- 
fective spacetime evolves dynamically as the BEC starts 
to expand. Then, the sonic horizon is formed as was 
shown by the above numerical simulation. If the effec- 
tive spacetime is static, we can introduce a following time 
coordinate: r = t + J vodr/(c 2 — Vq), and the effective 



spacetime metric becomes 



ds 2 oc -(c 2 - v 2 )dr 2 + 



l dr 2 



+ r 2 dn 2 s2 . (A2) 



From this expression, it is found that the horizon is lo- 
cated at the surface where the condition c s = \v \ is sat- 
isfied. A new coordinate v is introduced as v = r + r* 
where r* = J c s dr/(c 2 — Vq), which is a coordinate char- 
acterizing ingoing light-like (null) rays in the effective 
spacetime. 

We assume here that the initial state of the quantum 
field <p is the vacuum state for the static observer in the 
initial effective spacetime. Under the time evolution of 
the effective spacetime caused by the expansion of the 
condensate, the creation and annihilation operators for 
the field ip also evolve, and particle creation occurs. 

Now we consider an observer who moves along his 
or her outgoing geodesic with proper time A, crossing 
the horizon at A = 0. Hereafter, we term the observer 
geodesic observer. If we assume that the horizon is lo- 
cated at r — rji , the proper time A is related to the coor- 

dinate v there via A ~ — Aoe ~ v , where cjj = c s (rjj), 
a = 2cff<9 r (u — c s )| r=rH and Ao is a constant. The ingo- 
ing mode functions ip^ = e~ lu)v have A-dependence near 
the horizon as 



, - 2chUJ i i \\ 
exp i ln(— A) 

a 



(A3) 



Initially, the state is the vacuum for the static observer 
and therefore the geodesic observer would see no exci- 
tation at short distance, because there will be no much 
higher positive frequency excitations than those deter- 
mined by the time scale of the dynamical expansion of 
the BEC. If we ignore the short distance cut-off deter- 
mined by the healing length, or equivalently, if the lat- 
ter condition in Eq.([3|) is ignored, this A-dependence of 
the ingoing mode functions implies that the particle cre- 
ation from the horizon into the inside of the condensate 
has thermal spectrum with the temperature given by ([9]). 
Furthermore, even if the short distance cut-off is taken 
into account.it is known that the result does not change 
in principle [25|]. 

Therefore, the particle emission from horizon will oc- 
cur in the case of expanding condensate, even where the 
subsonic region is inside of the horizon. 
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